

#R version 4.0.2 
#Zelig version 5.1.7

##load packages
library(Zelig)
library(survival)
library(Amelia)
library(simPH)
library(car)

#load data (models 1 and 2)
load("ghstdata.RData")

#survival data (models 3 and 4)
load("ghstsurvdata.RData")

#survival tscs data Feb29 (models 5 and 6)
load("ghsttscssurvdatafeb29.RData")

#survival time to domestic measures (models 7 and 8)
load("ghstsurvdomesticmeasuresdata.RData")

#survival tscs data full period (Supp material Table S4.4 models )
load("ghsttscssurvdatafull.RData")

#Table 1

model1<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data=ghstdata, model="logit")
summary(model1)

model2<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data=ghstdata, model="logit")
summary(model2)

#Simulate first difference

set.seed(02474)
x.low<-setx(model2, natideol=quantile(na.omit(ghstdata$natideol))[2])
x.high<-setx(model2, natideol=quantile(na.omit(ghstdata$natideol))[4])
s.out1<-sim(model2,x=x.low, x1=x.high)
summary(s.out1)

##Table 2
###NOTE, standard errors for Cox models in article are of coefficient values not hazard rates exp(coef).

model3<- coxph(Surv(daysoutbreakstart,eventfeb29) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data = ghstsurvdata[which(ghstsurvdata$schengen!=1),])
summary(model3)

model4<- coxph(Surv(daysoutbreakstart,eventfeb29) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data = ghstsurvdata[which(ghstsurvdata$schengen!=1),])
summary(model4)

model5<- coxph(Surv(start,end,censor) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan,robust=TRUE, cluster=COUNTRY_CODE,
              data = ghsttscssurvdatafeb29[which(ghsttscssurvdatafeb29$schengen!=1),])
summary(model5)

model6<- coxph(Surv(start,end,censor) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion+casesbycountry+regionbarriers,robust=TRUE, cluster=COUNTRY_CODE,
              data = ghsttscssurvdatafeb29[which(ghsttscssurvdatafeb29$schengen!=1),])
summary(model6)

# Simulate hazard ratios
         
 set.seed(02474)        
 Sim1 <- coxsimLinear(model4, b = "natideol", Xj = quantile(na.omit(ghstsurvdata$natideol))[[4]], Xl=quantile(na.omit(ghstsurvdata$natideol))[[2]], nsim = 100, qi = "Hazard Ratio", ci=.95)
 
mean(Sim1$sims[3][[1]])
range(Sim1$sims[3])

set.seed(02474)        
 Sim1 <- coxsimLinear(model4, b = "distwuhan", Xj = quantile(na.omit(ghstsurvdata$distwuhan))[[2]], Xl=quantile(na.omit(ghstsurvdata$distwuhan))[[4]], nsim = 100, qi = "Hazard Ratio", ci=.95)
 
mean(Sim1$sims[3][[1]])
range(Sim1$sims[3])

#Table 3
###NOTE, standard errors for Cox models in article are of coefficient values not hazard rates exp(coef).

model7<- coxph(Surv(diffstayathomeandrestrict,censordiffstayathomeandrestrict) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data = ghstsurvdomesticmeasuresdata[which(ghstsurvdomesticmeasuresdata$schengen!=1),])
summary(model7)

model8<- coxph(Surv(diffstayathomeandrestrict,censordiffstayathomeandrestrict) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data = ghstsurvdomesticmeasuresdata[which(ghstsurvdomesticmeasuresdata$schengen!=1),])
summary(model8)

#Simualate hazard ratio

 set.seed(02474)        
 Sim1 <- coxsimLinear(model8, b = "natideol", Xj = quantile(na.omit(ghstsurvdomesticmeasuresdata$natideol))[[4]], Xl=quantile(na.omit(ghstsurvdomesticmeasuresdata$natideol))[[2]], nsim = 100, qi = "Hazard Ratio", ci=.95)
 
mean(Sim1$sims[3][[1]])

range(Sim1$sims[3])

#Figure 1

set.seed(02474)

seq<-seq(0, 1, .1)
ev.est=c()
ci.low=c()
ci.high=c()

for (i in seq){
  model2<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data=ghstdata, model="logit")
  x.low<-setx(model2, natideol=i)
  s.out1<-sim(model2,x=x.low)
  ev<-s.out1$get_qi(qi = "ev", xvalue = "x")
  ev.est<-cbind(ev.est, mean(ev))
  ci.low<-cbind(ci.low, quantile(ev, .025))
  ci.high<-cbind(ci.high,quantile(ev, .975))
  
}

par(mai=c(.7,0.7,.1,.2), mgp=c(2.5,1,0))
plot(seq,ev.est,xaxt='n', yaxs="i",yaxt='n', cex.main=0.70, cex.lab=1, main="" , lwd=1.5, type="n", ylim=c(0,1.05),xlim=c(0,1),xlab="", ylab="")
for (i in seq){
  lines(x=c(i,i),c(ci.low[seq==i], ci.high[seq==i]), type='l', lty=1, lwd=2.5, col="gray58" )}
points(seq,ev.est, pch= 19, cex=.8)
#abline(h=0,lty=2, lwd=.5)
axis(1,seq(0, 1, .1), cex.axis=.8)
axis(2,seq(-1,1, .10), cex.axis=.8)

mtext("Likelihood of imposing border restrictions", side=2, line=2.5, cex=.8)
mtext("Extent to which government is characterized by a nationalist ideology", side=1, line=2.5, cex=.8)

legend(.6, .3, c("95% confidence interval"), cex=0.8, lwd=c(2.5), lty=c(1), col=c("gray58"), bty="n")

#Supplementary material

#summary stats

ghstsurvdata$logtravtour<-log(ghstsurvdata$travtourshare)

summarydata1<-ghstsurvdata[c("restrictfeb29", "daysoutbreakstart","natideol", "polity", "ruleoflaw", "ghsiscore","loggdppc")]
summary(summarydata1)

summarydata1b<-ghstsurvdata[c("distwuhan", "whoregion", "gov_effect", "logpop", "air_travel", "logtravtour", "beltandroad", "globalizationindex")]
summary(summarydata1b)


###Table S4.1

#1. Border orientation

bo<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion+borderorientation, data=ghstdata, model="logit")
summary(bo)

#2. Past barriers

pastb<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion+h1n1orebolab, data=ghstdata, model="logit")
summary(pastb)

#3. Welcoming score
welcome<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion+welcomingscore, data=ghstdata, model="logit")
summary(welcome)

#4. Nationalist regime x democracy

interact<-zelig(restrictfeb29~natideol*polity+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data=ghstdata, model="logit")
summary(interact)

#5. Health capacity x democracy

demochi<-zelig(restrictfeb29~natideol+ghsiscore*polity+loggdppc+gov_effect+distwuhan+logpop+air_travel+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data=ghstdata, model="logit")
summary(demochi)

#4. Alternative measures

altm<-zelig(restrictfeb29~natideol+hdi+healthexpgdp+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+logtradegdp+logchinatradedep+beltandroad+whoregion, data=ghstdata, model="logit")
summary(altm)

#Table S4.2

cs<- coxph(Surv(daysoutbreakstart,event) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data = ghstsurvdata[which(ghstsurvdata$schengen!=1),])
summary(cs)

csfull<- coxph(Surv(daysoutbreakstart,event) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data = ghstsurvdata[which(ghstsurvdata$schengen!=1),]) 
summary(csfull)

ts<- coxph(Surv(start,end,censor) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan,robust=TRUE, cluster=COUNTRY_CODE,
              data = ghsttscssurvdatafull[which(ghsttscssurvdatafull$schengen!=1),])
summary(ts)

tsfull<- coxph(Surv(start,end,censor) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion+casesbycountry+regionbarriers,robust=TRUE, cluster=COUNTRY_CODE,
              data = ghsttscssurvdatafull[which(ghsttscssurvdatafull$schengen!=1),])
summary(tsfull)

#Table S4.3 Multiple imputation

ghstdata$logtravtourshare<-log(ghstdata$travtourshare)

data<-ghstdata[c("restrictfeb29", "natideol", "ghsiscore", "loggdppc", "gov_effect", "distwuhan", "logpop", "air_travel", "polity", "populist", "ruleoflaw", "logtravtourshare", "beltandroad", "globalizationindex", "whoregion")]

#set variable bounds 
bds <- matrix(c(9,-10,10,10, 0, 1, 13, 0, 1), nrow = 3, ncol = 3, byrow=TRUE)

#set seed
set.seed(02474)

#impute
a.out1<-amelia(data, m=5,idvars=c("whoregion"), ords=c("restrictfeb29", "polity", "populist", "beltandroad"),p2s= 1, bound=bds, empri=.01*nrow(data))

impu1<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data=a.out1, model="logit")
summary(impu1)

impu2<-zelig(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+logtravtourshare+beltandroad+globalizationindex+whoregion, data=a.out1, model="logit")
summary(impu2)

#Table S4.4

#Alternative DV for domestic measures models

adc1<- coxph(Surv(diffinternalrestrictandrestrict,censordiffinternalrestrictandrestrict) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data = ghstsurvdomesticmeasuresdata[which(ghstsurvdomesticmeasuresdata$schengen!=1),])
summary(adc1)

adc2<- coxph(Surv(diffinternalrestrictandrestrict,censordiffinternalrestrictandrestrict) ~ natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data = ghstsurvdomesticmeasuresdata[which(ghstsurvdomesticmeasuresdata$schengen!=1),])
summary(adc2)

#S4.5

##vif check

#Table 1 Models 1 and 2
vif1<-glm(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan, data=ghstdata, family="binomial")
summary(vif1)

vif2<-glm(restrictfeb29~natideol+ghsiscore+loggdppc+gov_effect+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+globalizationindex+whoregion, data=ghstdata, family="binomial")
summary(vif2)

vif(vif1)
vif(vif2)

#removing who region, gov effectiveness, globalization from Model 2 because of large VIF

vif3<-glm(restrictfeb29~natideol+ghsiscore+loggdppc+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad, data=ghstdata, family="binomial")
summary(vif3)

vif(vif3)


#removing same variables from other full models in main analysis to ensure results are consistent

#Table 2 model 4 and 6

vif4<- coxph(Surv(daysoutbreakstart,eventfeb29) ~ natideol+ghsiscore+loggdppc+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad, data = ghstsurvdata[which(ghstsurvdata$schengen!=1),])

summary(vif4)

vif5<- coxph(Surv(start,end,censor) ~ natideol+ghsiscore+loggdppc+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad+casesbycountry+regionbarriers,robust=TRUE, cluster=COUNTRY_CODE,
              data = ghsttscssurvdatafeb29[which(ghsttscssurvdatafeb29$schengen!=1),])
summary(vif5)

#Table 3 model 8

#domestic model

vif6<- coxph(Surv(diffstayathomeandrestrict,censordiffstayathomeandrestrict) ~ natideol+ghsiscore+loggdppc+distwuhan+logpop+air_travel+polity+populist+ruleoflaw+log(travtourshare)+beltandroad, data = ghstsurvdomesticmeasuresdata[which(ghstsurvdomesticmeasuresdata$schengen!=1),])
summary(vif6)


